Molecular hydrogen in graphite: A path-integral simulation 
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Molecular hydrogen in the bulk of graphite has been studied by path-integral molecular dynamics 
simulations. Finite-temperature properties of H2 molecules adsorbed between graphite layers were 
analyzed in the temperature range from 300 to 900 K. The interatomic interactions were modeled 
by a tight-binding potential fitted to density-functional calculations. In the lowest-energy position, 
an H2 molecule is found to be disposed parallel to the sheets plane. At finite temperatures, the 
molecule explores other orientations, but its rotation is partially hindered by the adjacent graphite 
layers. Vibrational frequencies were obtained from a linear-response approach, based on correlations 
of atom displacements. For the stretching vibration of the molecule, we find at 300 K a frequency 
uj s — 3916 cm' 1 , more than 100 cm -1 lower than the frequency corresponding to an isolated H2 
molecule. Isotope effects have been studied by considering also deuterium and tritium molecules. 
For D2 in graphite we obtained uj a = 2816 cm -1 , i.e., an isotopic ratio cj s (H)/oj s (D) = 1.39. 

PACS numbers: 61.72.S-, 63.20.Pw, 81.05.uf, 71.15.Pd 
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I. INTRODUCTION 

In the last few years there has been a surge of interest 
on carbon-based materials, specially on those composed 
of C atoms displaying sp 2 hybridization. This group of 
materials includes carbon nanotubes and fullerenes, dis- 
covered in last decades, as well as graphene, found in re- 
cent years' and the well-known graphite. These mate- 
rials, apart from its interest in basic science, are promis- 
ing tools for diverse technological applications. Thus, 
carbon-based systems, in general, are considered as pos- 
sible candidates for hydrogen storage^ Also, chemisorp- 
tion on two-dimensional systems, such as graphene or 
graphite surfaces, is supposed to be important for cat- 
alytic processes £ 

The scientific and technological interest of hydrogen 
as an impurity in solids and on surfaces has existed for 
several decades. In principle, it seems to be one of the 
simplest impurities, but a clear understanding of its prop- 
erties is not trivial because of its low mass, and needs 
the combination of advanced experimental and theoret- 
ical methodsi£i£ In addition to its basic interest as an 
impurity, a remarkable aspect of hydrogen in solids and 
surfaces is its capability of passivating defects and form- 
ing complexes, facts that have been extensively studied 
for many years 

Experimental investigations of isolated hydrogen in 
graphite turn out to be difficult because of the large sen- 
sitivity required to detect this impurity, along with the 
presence of a large amount of hydrogen trapped at the 
boundaries of graphite crystallites The stable config- 
urations of hydrogen in the bulk of graphite have been 
studied in several theoretical worksji 2 . - — where special 
emphasis was laid upon both atomic and molecular forms 
of this impurity. Moreover, theoretical techniques have 
been applied to investigate the diffusion, trapping, and 
recombination of hydrogen on a graphite surfaceJ^r— In 
this respect, chemisorption of a single hydrogen atom on 



a graphene sheet has been studied by several authors us- 
ing ab-initio methods^ii2r— and their results show the 
appearance of a defect-induced magnetic moment, along 
with a large structural distortion! 20 ! 22 ' 23 

For the storage of hydrogen in graphite one should also 
consider the presence of H2 molecules in the graphite 
bulk, which are expected to be physisorbed in the inter- 
layer space d 13 ' 14 Here we will focus on isolated hydrogen 
molecules trapped between graphite sheets. The impor- 
tance of this problem is twofold: it is interesting as a 
point defect in materials physics, for its relevance in the 
stability and diffusion of hydrogen in carbon-based solids, 
and also H2 in graphite is an example of a light molecule 
sitting and moving in a confined geometry, where quan- 
tum effects can be nontrivial. 

Earlier theoretical investigations of molecular hydro- 
gen in solids have focused on finding the lowest-energy 
position and stretching frequency of the molecule, in- 
cluding sometimes anharmonic effects obtained from the 
calculated potential-energy surface ) 24 " 28 and the quan- 
tum rotation of H2 molecules! 29 ' 30 Density-functional 
electronic-structure calculations in condensed matter are 
nowadays very reliable, but they usually deal with atomic 
nuclei as classical particles, so that typical quantum ef- 
fects like zero-point vibrations are not directly included. 
These effects can be taken into account by making use of 
harmonic or quasiharmonic approximations, but are dif- 
ficult to consider when large anharmonicities are present, 
as may happen for light impurities such as hydrogen. 

The quantum character of the atomic nuclei can be 
taken into account by using the path-integral molecular 
dynamics (or Monte Carlo) approach, which has been 
shown to be very useful in this respect. A notable bene- 
fit of this procedure is that all nuclear degrees of freedom 
can be quantized in an efficient and direct way, so that 
both quantum and thermal fluctuations are directly in- 
cluded in the calculations. Thus, molecular dynamics or 
Monte Carlo sampling applied to evaluate path integrals 
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allows one to carry out quantitative and nonperturbative 
studies of anharmonic effects in many-body systems! 31 ' 32 

In the present paper, we use the path-integral molec- 
ular dynamics (PIMD) method to study hydrogen 
molecules adsorbed in the interlayer region of graphite. 
Particular emphasis was placed upon anharmonic effects 
in their quantum dynamics at different temperatures. We 
analyze the isotopic effect on structural and vibrational 
properties of these molecules, by considering also molec- 
ular deuterium (D2) and tritium (T2). Path- integral 
methods similar to that employed in this work have been 
applied earlier to investigate hydrogen in metals 3 - 1 - and 
semiconductors j&M. as well as on surfaces! 38 ' 39 In rela- 
tion to the behavior of molecular hydrogen in confined 
regions, H2 has been studied inside carbon nanotubcs 
by diffusion Monte Carlo42 Also, path-integral simula- 
tion methods have been thoroughly applied to study con- 
densed phases of hydrogen in molecular form4i~— 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the models em- 
ployed in our calculations. Our results are presented 
in Sec. Ill, dealing with the spatial derealization of H 
atoms, interatomic distance, vibrational frequencies, and 
kinetic energy. Sec. IV includes a summary of the main 
results. 



II. COMPUTATIONAL METHOD 
A. Path-integral molecular dynamics 

Our calculations are based on the path-integral for- 
mulation of statistical mechanics. In this formulation, 
the partition function is evaluated by a discretization 
of the density matrix along cyclic paths, made up of a 
finite number L (Trotter number) of "imaginary-time" 
steps i 45 ' 46 In the implementation of this procedure to 
numerical simulations, such a discretization gives rise to 
the appearance of L "beads" for each quantum particle. 
These beads can be treated in the calculations as classi- 
cal particles, since the partition function of the original 
quantum system is isomorph to that of a classical one. 
This isomorphism is obtained by replacing each quantum 
particle by a ring polymer consisting of L classical parti- 
cles, connected by harmonic springs. 31 ' 32 In many-body 
problems, the configuration space can be adequately sam- 
pled by molecular dynamics or Monte Carlo techniques. 
Here, we have used the PIMD method, which was found 
to require less computer time for the present problem. 
We have employed effective algorithms for performing 
PIMD simulations in the canonical NVT ensemble, as 
those described in the literature ! 47 ' 48 

Calculations have been performed within the Born- 
Oppenheimer approximation, which allows us to define 
a 3iV-dimensional potential-energy surface for the mo- 
tion of the atomic nuclei. An important issue in this 
type of calculations is the proper description of inter- 
atomic interactions, which should be as realistic as pos- 



sible. Since effective classical potentials present many 
limitations to reproduce the many-body energy surface, 
one should resort to self-consistent quantum-mechanical 
methods. However, density functional (DF) or Hartree- 
Fock-based self-consistent potentials require computing 
resources that would appreciably restrict the size of our 
simulation cell and/or the number of simulation steps. 
We found a reasonable compromise by obtaining the 
Born-Oppenheimer surface from a tight-binding (TB) ef- 
fective Hamiltonian, derived from DF calculations ~ The 
capability of TB methods to simulate different properties 
of solids and molecules has been reviewed by Goringe et 
al£2. In particular, the ability of our DF-TB potential to 
predict frequencies of C-H vibrations in molecular sys- 
tems was shown in Refs. mini We have employed earlier 
this TB Hamiltonian to describe hydrogen-carbon inter- 
actions in diamon d 36 ' 37 and graphene.— The TB energy 
consists of two parts; one of them is the sum of energies of 
occupied one-electron states, and the other corresponds 
to a pairwise interatomic potential^ 9 - Since a reliable de- 
scription of the hydrogen molecule is essential for our 
purposes, particular attention was put on the H-H pair 
potential, which has been taken as in our earlier study 
of molecular hydrogen in silicon^ 3 - This pair potential re- 
produces the main features of known effective interatomic 
potentials for H2, such as the Morse potential^ 

Simulations were carried out on a graphite supercell 
containing 64 C atoms and one hydrogen molecule (H2, 
D2, or T2), and periodic boundary conditions were as- 
sumed. The simulation cell includes two graphite sheets, 
each one being a 4 x 4 graphene supercell of size 4a = 
9.84 A. An AB layer stacking was considered, so that 
both sheets are disposed in such a way that the center of 
each hexagonal ring of one of them lies over a C atom of 
the adjacent sheet. To keep this type of stacking along 
a simulation run, avoiding diffusion of the graphite lay- 
ers, the center-of-gravity of each layer was not allowed 
to move on the layer plane, which will be referred in the 
sequel as the (2, y) plane. The average distance between 
sheets is a half of the supercell parameter along the z axis 
(perpendicular to the graphite layers), and was taken to 
be 3.35 A. For the reciprocal-space sampling we have used 
only the T point (k = 0), since the effect of employing a 
larger k set is a nearly constant shift in the total energy, 
with little influence on the energy differences between 
different atomic configurations. The influence of the cell 
size on the results of the simulations has been checked 
by considering graphite supercells including up to 144 C 
atoms (a 6 x 6 supercell). The results found for 4x4, 
5x5, and 6x6 supercells coincided within statistical 
error bars. In particular, we checked the kinetic energy 
of the H2 molecule (error bar of ±3 meV) and the mean 
H-H distance (error bar of ±2 x 10 -4 A). Also, including 
more graphite layers in the simulation cell does not affect 
the results of the PIMD simulations. 

Sampling of the configuration space has been carried 
out at temperatures between 300 and 900 K. For com- 
parison, we also carried out PIMD simulations of pure 
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graphite, as well as simulations of hydrogen molecules 
between rigid graphite sheets (in which the C atoms are 
kept fixed on their unrelaxed positions; see Sect. II. B for 
a precise definition of this approach). Moreover, some 
simulations were performed in the classical limit, which 
is obtained in our context by setting the Trotter num- 
ber L = 1. The electronic-structure calculations were 
performed without considering a temperature-dependent 
Fermi filling of the electronic states, which is reasonable 
for the temperature range under consideration. For a 
given temperature, a typical simulation run consisted of 
2 x 10 4 PIMD steps for system equilibration, followed by 
10 6 steps for the calculation of ensemble average proper- 
ties. To keep roughly a constant precision in the PIMD 
results at different temperatures, the Trotter number was 
scaled with the inverse temperature (L cx 1/T), so that 
LT = 18000 K, which translates into L = 60 for T = 300 
K. Quantum exchange effects between hydrogen nuclei 
were not taken into account, as they are negligible at the 
temperatures considered here, and both atomic nuclei in 
a molecule were treated as distinguishable particles. 

The simulations were performed by using a staging 
transformation for the bead coordinates. The canonical 
ensemble was generated by coupling chains of four Nose- 
Hoover thermostats to each degree of freedom^ To inte- 
grate the equations of motion, we employed a reversible 
reference-system propagator algorithm (RESPA), which 
allows one to define different time steps for the integra- 
tion of fast and slow degrees of freedom.4 7 - The time 
step At associated to the calculation of DF-TB forces 
was taken in the range between 0.1 and 0.3 fs, which 
was found to be appropriate for the interactions, atomic 
masses, and temperatures considered here. For the evo- 
lution of the fast dynamical variables, associated to the 
thermostats and harmonic bead interactions, we used a 
smaller time step St = At/4. Note that for H2 in graphite 
at 300 K, a simulation run consisting of 10 6 PIMD steps 
needs the calculation of energy and forces with the TB 
code for 6 x 10 7 configurations, which required the use of 
parallel computing. 



along a PIMD simulation run is then given by: 

A ' = l(|>-< ? >) 2 )' ( 2 ) 

where (...) indicates a thermal average at temperature T. 
After a direct transformation, one can write A 2 as 

A, 2 = Ql + Cl , (3) 

with 
and 

C, 2 = ((r-(r)) 2 ) = (r 2 )-(r) 2 . (5) 

The first term, Q 2 , is the mean-square "radius-of- 
gyration" of the ring polymers associated to the quan- 
tum particle under consideration, and gives the average 
spatial extension of the paths£L The second term on the 
r.h.s. of Eq. (J3j) , C 2 , is the mean-square displacement 
of the path centroid. This term is the only one remain- 
ing in the high-temperature (classical) limit, since then 
each path collapses onto a single point and the radius- 
of-gyration vanishes. In cases where the anharmonicity 
is not very large, the spatial distribution of the centroid 
r is similar to that of a classical particle moving in the 
same potential, and C 2 can be considered as a kind of 
semiclassical derealization. 

As indicated above, results of our PIMD simulations 
for H2 in graphite will be compared with those obtained 
for rigid graphite layers. This means that in the latter 
case the centroids of the quantum paths corresponding to 
carbon atoms are kept fixed on their ideal atomic posi- 
tions, so that no relaxation of the host atoms is allowed in 
the presence of the hydrogen molecule. This restriction 
allows however paths of the C atoms to be delocalized 
around their ideal sites, i.e. in this approach one has 
C 2 = and Q 2 > for the carbon atoms. 



B. Path centroid derealization 

We now define some spatial properties of the particle 
paths that will be used in the analysis of the simulation 
results. The center-of-gravity (centroid) of the quantum 
paths of a given particle is calculated as 

? =zX>' (1) 

i=l 

Ti being the position of bead i in the associated ring 
polymer. 

The mean-square displacement of a quantum particle 



C. Anharmonic vibrational frequencies 

Vibrational frequencies are often employed as finger- 
prints of impurities in solids, revealing information on 
the position that they occupy and on their interactions 
with the nearby hosts atoms. A traditional approach for 
calculating vibrational frequencies of impurities consists 
in obtaining the eigenvalues of the dynamical matrix as- 
sociated to the atoms in the simulation cell, which yields 
the frequencies in a harmonic approximation. However, 
for light impurities the anharmonicity can be large, and 
the harmonic frequencies are only a first (maybe poor) 
approximation. 

Anharmonic frequencies of vibrational modes will be 
calculated here by using a method based on the linear 
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response (LR) of the system to vanishingly small forces 
applied on the atomic nuclei. In the context of path- 
integral simulations, this approach has been shown to 
represent a significant improvement with respect to a 
standard harmonic approximation^ In particular, the 
vibrational frequency of the H2 stretching mode is de- 
rived by the LR method as 



0.03 



(6) 



where kg is Boltzmann's constant, fi is the reduced mass 
of the H2 molecule, and C\ is the mean-square displace- 
ment of the H-H distance, d, that is obtained by a relation 
analogous to Eq. ([5]) , after substitution of the particle co- 
ordinate r by the interatomic distance d. Details on this 
method and discussions of its capability for predicting 
vibrational frequencies of molecules and solids are given 
elsewhere£i£i£&S 



III. RESULTS 
A. Atomic delocalization 

For an H2 molecule we find a minimum-energy position 
at an interstitial site between a carbon atom in a graphite 
sheet and an hexagonal ring in an adjacent sheet. At this 
position, the preferred orientation of the molecule is par- 
allel to the graphite planes, in agreement with earlier cal- 
culations based on density-functional theory— At finite 
temperatures the molecule will explore other positions 
and orientations with respect to the graphite layers. In 
particular, the molecule can be tilted, forming an angle 
ip with the (x, y) plane. 

In Fig. 1 we present the probability distribution of the 
angle if, as derived from our PIMD simulations at two 
temperatures: 300 K (solid line) and 900 K (dashed line). 
This distribution has a maximum at if = (H-H parallel 
to the layers) , and vanishes for H-H perpendicular to the 
sheet plane [if = ±90°). As temperature increases, the 
probability distribution broadens slightly, but it remains 
as a peak centered at if = 0. However, we find that the 
molecule is free to rotate in the (x, y) plane. 

Even though there is no direct bond between molecular 
hydrogen and the graphite layers, the later relax slightly 
in the presence of the H2 molecules, so that they follow 
the H2 motion in the interlayer space. This means that 
the molecules are more mobile in the presence of flexible 
layers than in the case of stiff graphite sheets, in which 
the C atoms are fixed on their unrelaxed (ideal) positions. 
This can be visualized by looking at the distribution of 
the angle ip in both cases at a given temperature. Thus, 
in Fig. 2 we display this probability distribution at T 
— 300 K. The dashed line corresponds to rigid graphite 
sheets, whereas the solid one was obtained for flexible 
sheets (mobile C atoms). As expected, the distribution 
of the angle ip is broader for flexible sheets, since in this 
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FIG. 1: Probability distribution of the angle ip between the 
H-H direction and the graphite sheets, as derived from PIMD 
simulations at two temperatures: 300 K (solid line) and 900 
K (dashed line). 
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FIG. 2: Probability distribution of the angle <p between the 
H-H direction and the graphite sheets, as derived from PIMD 
simulations at 300 K. Shown are results obtained for mobile C 
atoms (no restrictions, solid line), as well as from a simulation 
in which the C atoms where kept fixed on their unrelaxed 
positions (fixed centroids, dashed line). 



case the H2 molecule can adopt configurations that are 
inaccessible in the presence of rigid graphite layers. We 
have also calculated the same probability distribution for 
D2 and T2 for flexible sheets, and at 300 K it turns out 
to be similar to that shown in Fig. 2 for H2 (solid line), 
but slightly narrower. 

We now turn to study the spatial delocalization of each 
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FIG. 3: Spatial delocalization of atomic nuclei (protons) in 
H2 along the z coordinate (perpendicular to the sheet plane) , 
as a function of temperature. Circles represent the mean- 
square displacement of the centroid of the quantum paths, C z , 
and squares correspond to the mean-square radius-of-gyration 
of the paths, Q 2 Z . The total delocalization, A^, is shown by 
triangles. Error bars of the total and centroid delocalization 
are in the order of the symbol size, whereas those of the radius- 
of-gyration are less than the symbol size. Dashed lines are 
guides to the eye. 



atomic nucleus in hydrogen molecules, that is expected 
to include a nonnegligible quantum contribution. For our 
problem of H2 in graphite, we have calculated separately 
both terms giving the atomic delocalization in Eq. ([3]), 
for each atom in the molecule. For a given temperature, 
the term C 2 does not converge to a well-defined value 
along a PIMD simulation, due to the onset of molecular 
diffusion in the interlayer space. However, its component 
C 2 along the z axis is an equilibrium property of the 
molecule, as in fact it cannot diffuse across graphite lay- 
ers. In Fig. 3 we display the values of Q 2 (spreading of 
the quantum paths, squares) and C 2 (centroid delocaliza- 
tion, circles), as derived from our PIMD simulations for 
the H 2 molecule at several temperatures. The total spa- 
tial delocalization along the z coordinate, A|, is shown as 
triangles. In this plot, one observes that C 2 is larger than 
Q 2 in the whole temperature range considered. From the 
spatial delocalization Q 2 shown in Fig. 3, one can esti- 
mate an effective frequency for hydrogen motion in the 
z direction. In fact, in a harmonic approximation Q 2 
can be expressed analytically as a function of frequency 
and temperature , 31 i 58 and comparing the delocalization 
expected for different frequencies with that given by our 
PIMD simulations, we found an effective frequency uj z of 
about 1200 cm- 1 . 

For the spreading of the quantum paths of each H atom 
we obtain at room temperature Q 2 = 9.4 x 10~ 3 A 2 , and 
it decreases as temperature is raised. It is interesting to 



compare this value with that found for Q 2 in the case of 
unrelaxed graphite layers. When the C atoms are fixed 
on their ideal positions, we find at 300 K, Q 2 = 5.9 x 10~ 3 
A 2 , much lower than that found when the C atoms are 
allowed to relax in the presence of the hydrogen molecule. 
It is also interesting to compare these values of the quan- 
tum delocalization in the direction perpendicular to the 
graphite sheets, with that on the (x, y) plane. Our sim- 
ulations yield Q 2 = Q 2 = 1.07 x 10~ 2 A 2 and 9.9 x 10" 3 
A 2 , for free and fixed C atoms, respectively, indicating 
that the quantum motion of hydrogen on the (x, y) plane 
is only slightly affected by motion of the carbon atoms. 
In fact, diffusion of the H2 molecules in the interlayer 
space occurs basically by classical jumps, as described 
elsewhere.— However, the quantum motion in the z di- 
rection is markedly affected by the relaxation of the C 
atoms, that contributes to enhance Q 2 by a factor of 1.6. 
In other words, this increase in the quantum delocaliza- 
tion is associated to a decrease in frequency (softening) of 
the vibrational modes of the H 2 molecule when full mo- 
tion of the C atoms is taken into account. These modes 
correspond to a displacement of the whole molecule along 
the z axis, and to the frustrated molecular rotation, with 
changes in the angle cp shown in Fig. 1. We note that the 
quantum paths have an average extension of ~ 0.1 A at 
300 K, much smaller than the H-H distance, thus justi- 
fying the neglect of quantum exchange between protons. 

For the D 2 molecule we obtain at 300 K, Q 2 — 5.4 x 
10 -3 A 2 in the case of free motion of all atoms in the 
simulation cell. Comparing with the H 2 molecule, we 
have Q 2 (H)/Q 2 (D) = 1.7, clearly higher than the low- 
temperature limit in a harmonic approximation, given by 
a ratio of y/2. Note that in the high-temperature limit 
Q 2 goes to zero, but the ratio Q 2 (H)/Q 2 (D) converges 
to the inverse mass ratio ; 31 i 58 in this case ra^/mn — 
2. For T 2 we found at 300 K, Q 2 = 3.8 x 10~ 3 A 2 , 
so that Q 2 (H)/Q 2 (T) = 2.5, also between a ratio of \/3 
expected at low temperature in a harmonic approach, 
and the high-temperature limit given by tot/tih = 3. 



B. Interatomic distance 

We first present results for classical calculations at zero 
temperature, where the atoms are treated as point-like 
particles without spatial delocalization. The interatomic 
potential employed here gives reliable results for molecu- 
lar hydrogen in vacuo (an isolated H 2 molecule). In par- 
ticular, the lowest-energy molecular configuration corre- 
sponds to a distance i?o between hydrogen atoms of 0.741 
A. At this distance we obtain for H 2 in a harmonic ap- 
proximation a stretching frequency of 4397 cm . 

The interatomic distance between hydrogen atoms in- 
creases when the molecule is introduced from the gas 
phase into the graphite bulk, due to an attractive in- 
teraction between H and the nearby C atoms. For 
the minimum-energy distance we found in this case Rq 
= 0.753 A, which is similar to that found for the H 2 
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FIG. 4: Mean distance between atoms in a hydrogen 
molecule in graphite, as derived from PIMD simulations for 
Eh (squares), D2 (circles), and T2 (triangles). Results ob- 
tained from classical molecular dynamics simulations are also 
shown for comparison (diamonds). Error bars are in the or- 
der of the symbol size. Dashed lines are linear fits to the data 
points. 
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FIG. 5: Mean interatomic distance for H2 molecules in 
graphite, as a function of temperature. Squares indicate re- 
sults derived from PIMD simulations in which all H and C 
atoms are mobile, whereas circles correspond to simulations 
in which all C atoms where kept fixed on their unrelaxed po- 
sitions (fixed centroids). Dashed lines are linear fits to the 
data points. 



molecule in the bulk of semiconductor materials with the 
same interatomic potential^ This interatomic distance 
is expected to rise for increasing temperature. In fact, 
in a classical approximation we obtained R = 0.756 A 
and 0.761 A at 300 K and 900 K, respectively. These 
classical results are presented in Fig. 4 as diamonds, 
and display a linear temperature dependence with slope 
dR/dT = 8.8 x 10" 6 A/K. 

PIMD simulations can be also employed to study the 
temperature dependence of the mean interatomic dis- 
tance R in a quantum model. The molecular expansion 
with respect to the lowest-energy classical geometry is 
due to a combination of anharmonicity in the stretching 
vibration of the H2 molecule and a centrifugal contribu- 
tion caused by molecular rotation. At 300 K we find for 
H2 a mean interatomic distance R = 0.777 A, to be com- 
pared with R = 0.756 A in the classical limit at the same 
temperature. This means that the interatomic distance 
of H2 in graphite increases by 0.02 A (about a 3%) when 
quantum effects are considered. For the molecules D2 
and T2, one expects smaller interatomic distances due 
to their larger mass and smaller vibrational amplitudes. 
In fact, at 300 K we found for D 2 , R = 0.770 A (i.e., 
7 x 10" 3 A less than for H 2 ), and for T 2 , R = 0.767 A. 
In Fig. 4 we present the temperature dependence of the 
mean distance for H2 (squares), D2 (circles), and T2 (tri- 
angles), as derived from our PIMD simulations for full 
quantum motion of molecular hydrogen and host atoms. 
For D 2 and T 2 we find a slope dR/dT = 3.1 x 10" 6 
A/K and 2.8 x 10 -6 A/K, respectively, close to the value 
obtained for H 2 : dR/dT = 2.6 x 10" 6 A/K. Note that 



these changes of the interatomic distance derived from 
the PIMD simulations are much smaller than that found 
in the classical limit: dR/dT = 8.8 x 10~ 6 A/K. 

It is interesting to compare these changes in the mean 
distance R with those corresponding to molecular hydro- 
gen in the gas phase. With this purpose we carried out 
some PIMD simulations of an isolated hydrogen molecule 
with the same interatomic potential at several temper- 
atures. These simulations yielded an increase in R as 
temperature is raised, given by dR/dT = 7.5 x 10~ 6 
A/K, a value three times larger than that found for H2 
in graphite. 

It is also interesting to analyze the effect of the motion 
of the graphite layers on the H-H distance. We have 
calculated this distance in PIMD simulations in which 
the C atoms are kept fixed on their unrelaxed sites. The 
results of the interatomic distance are shown in Fig. 5 
as circles. We find that the H-H distance is in this case 
smaller than that obtained for full quantum motion of all 
the atoms in the simulation cell. In fact, at 300 K the 
distance R decreases by about 4 x 10 -3 A. However, the 
slope dR/dT for the fixed-lattice model is larger than in 
the case of mobile C atoms. In fact, for fixed C atoms 
we find dR/dT = 5.3 x 10" 6 A/K, to be compared with 
a slope of 2.6 x 10 -6 A/K obtained for flexible graphite 
sheets. 

From the zero-temperature classical calculations, we 
found that the H2 molecule is expanded when it is in- 
troduced into the graphite bulk, as a consequence of the 
attractive C-H interaction. At finite temperatures, the 
H-H distance is also controlled by the centrifugal effect 
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caused by molecule rotation. In fact, a 3D rotation is par- 
tially frustrated in the interlayer region, as shown above. 
In this respect, molecular rotation is favored by the re- 
laxation of graphite sheets, thus yielding a larger cen- 
trifugal contribution to the molecule expansion than for 
rigid sheets. 

The average interatomic distance allows us to estimate 
a moment of inertia for the molecule, and then the wave- 
number difference a>oi between J = and J = 1 ro- 
tational levels. This gives for H2, woi ~ HO cm -1 , 
somewhat lower than that known for the free molecule 
in vacuo— (wqi = 118.6 cm -1 ). For H 2 in graphite, 
however, the J = 1 level will be split due to the hin- 
dered motion of the molecule for H-H perpendicular to 
the graphite sheets. Unfortunately, the magnitude of this 
splitting is not accessible from our PIMD simulations. 



C. Stretching frequency 

The stretching frequency of H 2 is an important fin- 
gerprint of the molecule, that can be used to detect and 
characterize this impurity in solids. We first note that the 
lowest-energy configuration of an isolated H2 molecule al- 
lows us to predict a stretching frequency of 4397 cm -1 
in a harmonic approximation. From path-integral simu- 
lations combined with the LR approach presented above 
in Sect. II. C, we obtain for a single H2 molecule at 300 
K a frequency w s — 4055 ± 5 cm -1 , i.e., the anharmonic 
shift amounts to more than 300 cm -1 . Note that in the 
vibrational frequencies derived from the LR procedure, 
the error bars are due to the statistical uncertainty asso- 
ciated to the PIMD simulations. 

The stretching frequency uj s is further reduced when 
the H2 molecule is inserted between the graphite sheets. 
In fact, from our PIMD simulations at 300 K we found 
for H2 a stretching frequency ui s — 3916 ± 8 cm -1 . This 
frequency is found to decrease slightly as the temperature 
is raised, as shown in Fig. 6 (squares). A linear fit to the 
data points gives a slope duj s /dT = —0.03 cm -1 /K. 

The quantum treatment of atomic nuclei in molecu- 
lar dynamics simulations is decisive to give a reliable de- 
scription of the vibrational frequencies of light atoms like 
hydrogen. In fact, we have applied the LR method to cal- 
culate the stretching frequency ui s from classical simula- 
tions. At 300 K we found for H2 in graphite a frequency 
uj s = 4082 ± 7 cm" 1 (for full motion of interstitial hydro- 
gen and host atoms), about 160 cm" 1 higher than that 
found in the full quantum simulations (w s = 3916 cm" 1 ). 
In Fig. 6 we have plotted the frequency ui s derived from 
the classical molecular dynamics simulations in the tem- 
perature range from 300 to 900 K. The overestimation of 
vibrational frequencies in a classical approach, in com- 
parison with the quantum results is usual in this kind 
of simulations^ since the classical calculations tend to 
yield results much closer to the harmonic approximation, 
which gives in general frequencies higher than the anhar- 
monic ones (as is the case here). 
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FIG. 6: Frequency of the stretching vibration of the H2 
molecule in graphite as a function of temperature. Symbols 
represent results derived from PIMD simulations in two dif- 
ferent approaches: squares, full quantum motion of H and C 
atoms; circles, classical motion of H and C atoms. Dashed 
lines are linear fits to the data points. Error bars correspond 
to the statistical uncertainty in the molecular dynamics sim- 
ulations. 



As presented above when discussing the interatomic 
distance in the hydrogen molecule, there are two 
main factors controlling the stretching frequency of the 
molecule in the graphite bulk. The first one is the in- 
teraction with the graphite sheets, which tends to en- 
large the H-H distance, with a concomitant decrease in 
the frequency uj s . This is the main factor contributing 
to the reduction observed when comparing results of an 
isolated molecule and a molecule in the interlayer region. 
The second factor is the molecular rotation, which is par- 
tially hindered between the graphite sheets, but in gen- 
eral causes a decrease in the vibrational frequency due 
to ro-vibrational coupling. All together, we find that 
oj s decreases only slightly as T is raised. This contrasts 
with the results obtained for the stretching frequency of 
H2 in the interstitial space of silicon from PIMD simula- 
tions similar to those presented here^ In that case, the 
molecule is free to rotate in the silicon bulk, and uj s is 
found to decrease about eight times faster than for H 2 in 
graphite. 

For the D2 molecule in graphite we find at 300 K a 
stretching frequency u> s = 2816 ± 5 cm" 1 . For increasing 
temperature, we obtained a trend similar to that found 
for H2, with a linear decrease in w s . For the isotopic shift 
we found a rather constant ratio between the stretching 
frequencies of H 2 and D 2 , that amounts to 1.39, slightly 
smaller than the ratio expected in a harmonic approx- 
imation: cl> s (H)/oj s (D) = 1.41. Experimentally, a ratio 
of 1.39 has been observed for the frequencies of these 
molecules in the gas phase. 60 For comparison, we men- 
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FIG. 7: Temperature dependence of the kinetic energy of 
molecular hydrogen in graphite, as derived from PIMD sim- 
ulations. Open symbols indicate results derived from sim- 
ulations with full quantum motion of all atoms in the cell: 
squares for H2, circles for D2, and triangles for T2. For com- 
parison we also present results for H2 with fixed C atoms 
(filled squares). Error bars are on the order of the symbol 
size. Dashed lines are guides to the eye. The dotted line 
corresponds to the classical limit with six degrees of freedom 
(Ef = 3k B T). 



tion that in a classical simulation of D2 at 300 K we found 
a stretching frequency u> s — 2896 ± 5 cm -1 , which yields 
an isotopic ratio of 1.41, as in a harmonic approach. From 
PIMD simulations of T2 in graphite we obtained at 300 
K, w s = 2324 ± 5 cm" 1 , so that w,(H)/w a (T) = 1.69, 
slightly lower than the harmonic expectancy of 1.73. 



D. Kinetic energy 

Path integral simulations allow one to obtain the ki- 
netic energy Ek of the quantum particles under consid- 
eration, which is basically related to the spread of the 
quantum paths. In fact, for a particle at a given tem- 
perature, the larger the mean-square radius-of-gyration 
of the paths, Q^, the smaller the kinetic energy. Here we 
have calculated Ek by using the so-called virial estima- 
tor, which has an associated statistical uncertainty lower 
than the potential energy of the system i 55 i 61 

To analyze the kinetic energy associated to the defect 
complex, we calculate Ek for the simulation cell with and 
without the hydrogen molecule: Ek (defect) = -Efc(64 C 
+ H 2 ) - Ek(64 C), where we use results obtained in both 
series of PIMD simulations, with and without the hydro- 
gen molecule in the graphite cell. In Fig. 7 we display 
the kinetic energy as a function of temperature for H2 
(squares), D 2 (circles), and T 2 (triangles). As expected, 
Ek increases as temperature rises, and at a given tem- 



perature, it is larger for smaller isotopic mass. At 300 K, 
we find E k = 0.238 eV for H 2 , 0.181 eV for D 2 , and 0.154 
eV for T 2 , which gives isotopic ratios: Ek(H-2) / Ek(T>2) = 
1.31 and £?fc(rI 2 )/2?fc(T 2 ) = 1.55. These ratios decrease 
as temperature is raised, and at 900 K they amount to 
1.16 and 1.24, respectively, as in the high-temperature 
(classical) limit they should converge to unity. For com- 
parison we also present in Fig. 7 the kinetic energy cor- 
responding to a classical model with six degrees of free- 
dom [E^ = 3fcsT, dotted line). For rising tempera- 
ture, the classical kinetic energy approaches the results 
of PIMD simulations, in particular those corresponding 
to the heaviest isotope (tritium), but at 900 K it is still 
lower than i?fc(T 2 ) by 43 meV (about 15% of the quan- 
tum value). 

In the low-temperature limit the quotient 
Ek (H2)/-Efe (D2) is expected to be close to 1.41, as 
obtained in a harmonic approximation. From our PIMD 
simulations at 300 K we obtained a ratio clearly lower 
than this value, what can indeed be due to the presence 
of anharmonicities in the molecular motion, but more 
importantly to the excitation of quantum levels higher 
than the ground state at this finite temperature. This is 
particularly true for molecular rotation, which is found 
to be rather free in the (x, y) plane. Something similar 
can be said from the ratio -E/£(H 2 )/£^(T 2 ) at 300 K, 
which is clearly lower than 

We have also calculated the kinetic energy of the hy- 
drogen molecule between rigid graphite sheets. The re- 
sults for Ek obtained in this case are shown in Fig. 7 
as filled squares. At 300 K we find Ek — 0.361 eV, 
clearly higher than in the case of mobile carbon atoms. 
This decrease in kinetic energy of H 2 for flexible graphite 
sheets is mainly due to a softening of the vibrational 
modes corresponding to motion of the center-of-mass of 
the molecule in the interlayer space. In fact, relaxation 
of the C atoms in the presence of the H 2 molecule causes 
a decrease in the energy barriers confining the molecule 
at a given position, eventually favoring its motion in the 
graphite bulk. This decrease in Ek is then related to 
the larger vibrational amplitude of the whole molecule 
between flexible graphite layers, indicating a nonnegligi- 
ble coupling in the motion of interstitial molecule and 
host atoms. In connection with this, in an earlier work 
based on classical molecular dynamics simulations^ it 
was shown that relaxation of the graphite sheets in the 
presence of the H 2 molecule helps to enhance molecular 
diffusion in the interlayer space. In particular, it was 
found a behavior that could not be simply explained by 
a single activation energy, but suggested the presence of 
correlations between successive molecular hops. 



IV. SUMMARY 

We have presented results of PIMD simulations for iso- 
lated hydrogen molecules adsorbed in the interlayer re- 
gion of graphite. This kind of simulations allow us to 
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calculate kinetic and potential energies at finite temper- 
atures, taking into account the quantization of host-atom 
motions, which is not easy to consider in fixed-lattice cal- 
culations. This includes consideration of zero-point mo- 
tion of guest and host atoms, which can be coupled in a 
non-trivial way in the many-body problem. Also, isotope 
effects can be readily explored, since the impurity mass 
appears as a parameter in the calculations. 

Hydrogen molecules are found to be disposed basically 
parallel to the graphite-layer plane, and free to rotate 
in this plane. Although thermal and quantum dereal- 
ization allow the molecule to explore other orientations, 
molecular rotation is restricted by the nearest graphite 
sheets and in fact the H-H axis is not found to approach 
the direction perpendicular to the layers. 

An important feature of H2 molecules adsorbed in 
solids is their stretching vibration lj s . For molecular 
hydrogen in graphite, we find at 300 K a frequency uj s 
= 3916 cm" 1 , to be compared with lu s = 4055 cm -1 
obtained for an isolated molecule with the same inter- 
atomic potential and at the same temperature. For D2 in 
graphite we find a frequency of 2816 cm -1 , which gives 
an isotopic ratio u; s (H)/w s (D) = 1.39, similar to those 
measured for free hydrogen molecules. It is remarkable 
that classical simulations yield for H2 a frequency uj s 
about 160 cm -1 larger than the PIMD simulations. The 
stretching frequency of H2 and D2 is found to decrease 
slightly as temperature rises, as a consequence of cou- 
pling with molecular rotation and anharmonicities in the 
interatomic potential. 



Results of the PIMD simulations including full quan- 
tum motion of all atoms have been compared with those 
obtained for rigid graphite layers. This comparison has 
shown that motion of carbon atoms affects appreciably 
several properties of the adsorbed molecules, i.e., inter- 
atomic distance, stretching frequency, kinetic energy, and 
atomic derealization. 

A challenging question, that should be taken into ac- 
count in future work, refers to considering coupling be- 
tween nuclear spins in the hydrogen molecule, i.e. dealing 
separately with ortho and para-El^. This is particularly 
important at low temperatures, where the quantum na- 
ture of molecular rotation has to be explicitly consid- 
ered in the simulations. A quantum treatment of the full 
problem is not trivial, being mainly complicated by the 
ro-vibrational coupling. Apart from equilibrium PIMD 
simulations such as those presented here, one can ap- 
ply similar methods to study quantum diffusion of H2 in 
graphite, by calculating free-energy barriers as in the case 
of atomic hydrogen in metals^ and semiconductors! 3 ^ 62 
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